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We consider systems of a small number of interacting bosons confined to harmonic potentials 
in one and two dimensions. By exact numerical diagonalization of the many-body Hamiltonian we 
O , determine the low lying excitation energies and the ground state energy and density profile. We 

discuss the dependence of these quantities on both interaction strength g and particle number N . 
The ground state properties are compared to the predictions of the Gross-Pitaevskii equation, and 
the agreement is surprisingly good even for relatively low particle numbers. We also calculate the 
specific heat based on the obtained energy spectra. 
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I. INTRODUCTION 



Recent observations of Bose-Einstein condensation in alkali gases confined to harmonic potentials has led to an 
extensive study of the ground state and low energy properties of trapped interacting Bose gases. A major part of this 
work, see e.g. [K*^, presents calculations based on the Bogoliubov approximation or Hartree-Fock theory. Predicted 
' excitation frequencies agree well with experiments ]lO|JTl[ | . The energy spectrum has also been described in terms of a 
hydrodynamic formalism [p^ , which agrees with Bogoliubov theory at low energy |^Jl^. In the Bogoliubov approx- 
I imation the system is described in terms of a dominating condensate. Fluctuations induced by finite temperature 
and interactions are assumed small. At zero temperature, when the effects of excited particles are neglected, the 
Ph condensate satisfies the Gross-Pitaevskii equation 1 . 

^ : The Bogoliubov approximation is applicable for systems with a large number of particles. In this article we will, 
I I ' on the contrary, consider systems with relatively few particles, i.e. iV 10 — 40. Finite N effects are then important, 
\ and one should in principle use the complete machinery of many-body theory when describing these systems. Here, 
' we present the results of such calculations, limiting ourselves to one and two dimensions. The use of such time- and 
computer memory consuming methods puts severe limitations on the size of the system one may describe. On the 
other hand, it is for such systems the corrections to mean-field calculations are important. 

The results presented here are therefore not immediately relevant for the experiments performed with trapped Bose 
J — , gases so far. Our aim is rather to gain some insight into the effects of the approximations which form the basis of the 
' mean-field theories mentioned above, in the limit of small particle numbers. We are also interested in effects which 
C ■ , are specific to one and two dimensions. In addition, future experiments in lower dimensions may explore the range of 
' parameters used here. 

The plan of the article is as follows. In the next section we recall the many-body formalism and describe the steps 
leading to the diagonalization of the many-body Hamiltonian. The Gross-Pitaevskii approximation, to which we 



compare our results for the ground state, is discussed in Section III. Nimrerical results are presented and discussed 



in Section [V. Finally we draw some conclusions in the last section. The appendix discusses effective couplings in 
^ ' highly anisotropic systems. 

o ■ 



II. MANY-BODY DESCRIPTION 

As a model for a trapped interacting Bose gas we study the Hamiltonian density 

V2 



„ +K.(r) 
2 m 



*(r)-)-W,(*t,^,), (2.1) 



where ^' is a bosonic field operator and Vex is the external trapping potential. Units where Ti — 1 are used. We will 
only consider isotropic harmonic potentials of the form Vexi^) — ^muPr^. Hi describes the interaction. Assuming 
short-range two-body interaction, this may be written 

Tl!/ = .g*^^t^*(r). (2.2) 
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In the s-wave approximation the coupHng constant, or interaction strength g is related to the scattering length a 
by the equation g = 2iTa/m in three dimensions. In Appendix ^ we discuss modifications of this relation in highly 
anisotropic harmonic traps. 

Harmonic modes. As a step towards finding the energy eigenvalues, we expand the field operator in a complete 
set of modes 

*(r) =^(pfc(r)afc. (2.3) 

k 

Here annihilates a particle in the state k and a\, creates a particle in the same state. The operators satisfy 
the commutation relation [afc,a|] = 5k,i, with other commutators vanishing. In translation invariant systems it 
is convenient to choose a set of plane waves. Here we choose to use the harmonic oscillator cigenfunctions. The 
Hamiltonian then takes the form 

H = j drH ^ LUkalak + g ^ Iki7nna\a\ajnan. (2.4) 

k k,Lrn,n 

Here, u>k is the single particle energy of level k, and fkimn is an overlap integral over four oscillator cigenfunctions 

fklmn= J drifiliplipm'Pni'r)- (2.5) 

Notice that fkimn is symmetric under permutation of the indices. For the complete set of cigenfunctions {^Pk} we 
will choose Hermite or Laguerre polynomials for one and two dimensions, respectively. The corresponding overlap 
integrals can be calculated numerically, though some of them are tabulated JTzt . The interaction conserves parity. 
In one dimension this leads to the constraint k + l + m + n — even integer for non-zero /. In two dimensions the 
rotational symmetry of the external potential suggests the use of radial and angular quantum numbers. The inter- 
action conserves angular momentum, and this again puts a constraint on the possible combinations of creation and 
annihilation operators. 



Many particle basis and diagonalization. We apply the configuration interaction approximation and diagonalizc the 
Hamiltonian in a many particle basis, which in the occupation- number representation may be written in the form 

IV'q) = l'^0'^l«2 • ■ ■"-k-)q, (2.6) 

with a labeling the different distributions of particles. Normalization to unity is assumed. Here Uk is the particle 
number in the single particle state k. The particle number is conserved: — N . In a practical calculation we 

must truncate the basis set at some upper state, here denoted by K . This sets an upper limit for the single particle 
energy of the states considered. Such a truncation is however somewhat unnatural as far as the many particle energy is 
concerned, since it includes the many particle state |Oo • • • Nk) with energy Nujk, but not the state |(-/V— l)o ■ • ■ ^k+i) 
which only has energy ujk+i + [N —l)u!o. A more consistent way of truncating the basis set is therefore to include 
only those states which have a total energy up to some maximal value Emax- In the actual calculations, Emax will be 
raised in steps until convergence is reached. The creation operator and annihilation operator act as follows 

flfcl^o • • • nii) = y/n^luQ ■■ - Hk-l - ■■ uk), 

al\no - ■ ■ UK) ^ Vnk + l|no ■ ■ ■ Uk + I ■ ■ ■ uk) ■ (2.7) 

The number operator is thus a^a^ with eigenvalue n^. 

With a basis set at hand, we may calculate the Hamiltonian matrix elements 

Ha', a = {iJa'\H\i>oi) = HoaSa' .a q- (2.8) 

The free part of the Hamiltonian contributes to the diagonal terms with an amount w^nfe, while the interaction 
part yields the contribution 

IV-o). (2.9) 
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The symmetries of the interaction can be used to reduce the Hamiltonian matrix to a block diagonal form. Splitting 
the complete Fock space into several symmetry invariant subspaces reduces the dimensions of the Hamiltonian matrix 
considerably and enables numerical diagonalization of much larger systems. In one dimension there is an even and an 
odd parity subspace and the ground state is a superposition of even parity states. In two dimensions the subspaces are 
labeled by the total angular momentum of the many particle state and the ground state has zero angular momentum. 

The Hamiltonian matrix is relatively sparse and the iterative Lanczos method is thus well suited for finding the 
ground state and the lowest excited states. This method enables us to diagonalizc much larger matrices than would 
be possible with standard library routines. However, we use the latter when finding the complete set of energy levels 
needed for calculating thermodynamic quantities. 

We would here like to emphasize that we have checked the programs thoroughly: They are seen to reproduce well 
known results in the case of zero coupling. We have also done some calculations for smaller systems by hand, and 
find that both the generation of the Hamiltonian matrix and the diagonalization works correctly. 

Density profile. Once the Hamiltonian is diagonalized, we have the coefficients for the low lying states, and es- 
pecially for the ground state. This enables us to find the ground state density distribution p{v). The density operator 
may be written 

= ^|fc)nfc,(/|, (2.10) 

where \k) is a single particle state, and Uki — {ijj\a\a]J\'i})) . The density is then p(r) = (r|f;,^|r). For a state 1-0) = 
Z^a ^alV'o) this reads 

p(r) - ^^,>fc(r) C*M^Po.'\a\ak\i^^). (2.11) 

k^l a, a' 



Dipole mode. The excitations in this system will in general have interaction dependent energies. However, as discussed 
by Fetter and Rokhsar p3[ |, there exists a dipole mode for each spatial dimension of the trap, which corresponds to a 
harmonic oscillation of the center of mass of the condensate. In jl^ the corresponding raising operator was constructed 
in the first quantized formalism: 

N 

4,^T.blv (2-12) 

i=l 

where bj^^ is the raising operator for particle number i in dimension /3. (We use a different notation than in [ p^ , in 
order to avoid confusion with quantities defined here.) With a two-body potential of the form V{ri — rj) it is easily 
shown that the excitation energy is uj, i.e. independent of the interaction. Thus, for each energy level there must be 
a ladder of levels with energy spacing u above. These dipole modes have been found in calculations based on the 
Bogoliubov approximation and in the hydrodynamic description JT^ . The occurrence of these states will serve 
as a check of convergence of the calculated energy levels. For completeness we note that the raising operator in terms 
of creation and annihilation operators in one dimension takes the form 

=J2 Vfc + l4+iafc. (2.13) 

k 

Here k labels the Hermite polynomials. Commuting the raising operator with the free particle Hamiltonian we find 

[Ho,A^]^LoA\ (2.14) 

and the excitation energy is thus uj. The vanishing of the commutator of A''^ with the interaction term is less obvious 
in this formalism. The commutator becomes 

[Hi, A^] 25 ^ al^/l + lifkmnl+lClnai - fk (2.15) 

and for each combination a^a|ama„ the prefactor vanishes. This is due to the following relation between the overlap 
integrals: 



3 



Vd+ Ifabcd+l + Vc + Ifabc+ld — \fa,fa-lbcd + ^fab-lcd- 



(2.16) 



Breathing mode. In two dimensions the many-body Hamiltonian has an SO{2, l)-symmetry discovered by Pitaevskii 
and Rosch This symmetry gives rise to excitations of energy 2uj identified with the breathing mode of the 

condensate. The excitation spectra for the planar systems must therefore contain ladders both with energy spacing 
Lu and 2u!. The symmetry is broken when a UV energy cut-off is introduced in the system. Such a cut-off is inherent 
in our calculation. We therefore expect that the level spacing 2 lo will be less accurate than the spacing uj caused by 
the dipole excitations. 



III. MEAN FIELD APPROXIMATION 



One of the purposes of this paper is to compare the predictions of the many-body calculation with those of mean 
field theory for quantities such as the ground state energy and density distribution. The condensate is in the Gross- 
Pitaevskii (GP) approximation [ p^JTs) described by a classical, macroscopic Bose field governed by the Hamiltonian 
density per particle 

Hci/N = <P*ir) (-^ + VeAr)) $(r) + .giV ($*$(r))^ (3.1) 
\ 2m / 

The field is here rescaled to satisfy J (ir$*$ — 1, thus the factor N in front of the interaction term. Notice that the 
coupling g and the particle number N occur only in the combination gN, as opposed to the many-body description. 
Minimalization of the corresponding Hamiltonian Hd = J dvTLci leads to the non-linear Schrodinger (or Gross- 
Pitaevskii) equation 

- + Kx(r) - $(!■) + 2g7V$*$$(r) = 0. (3.2) 
Im J 

The chemical potential introduced is given by the normalization of 4>. This equation can be derived from the Hartree- 
Fock-Bogoliubov approximation when the effect of excited particles is neglected . 

Oscillator basis. The equation for $ may be solved introducing a harmonic oscillator basis 

<&(r) =^Cfc0fc(r). (3.3) 

k 

Normalization now corresponds to jc^p = 1. The equation then takes the form 

^(wfc - ^)cfc(/)/c(r) + 2gN ^ cldCmCplfjiicpniir) = 0. (3.4) 

fc k,Lm 

Multiplying with (^* and integrating we finally end up with the set of equations 

{uJn - n)Cn + 2gN ^ clciCmfklmn = 0, (3.5) 
fc,/,m 

where / again is the overlap integral over four harmonic oscillator eigenfunctions. The basis set is truncated at some 
level K, and K is raised until convergence is reached. The corresponding finite set of equations may be solved by 
use of the Newton-Raphson method. This method was applied in |0] for the case of a three-dimensional Bose gas. 
Solutions in one dimension have also been obtained earlier |^ . Again the symmetry of the ground state may be used 
to choose the proper subset of eigenfunctions. In one dimension we only have to consider even functions. In two 
dimensions only zero angular momentum eigenfunctions contribute. 

Once the coefficients are determined one must check if the normalization condition is satisfied. If not, a dif- 
ferent chemical potential is chosen, and the process is repeated until normalization is obtained. Having the correct 
coefficients, we easily find the ground state density profile 
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p{v) = \Y,CkM^)\'. (3.6) 

k 

In addition the ground state energy per particle reads 

k klmn 



IV. NUMERICAL RESULTS 



In the following we present the results of the numerical diagonalization of the many-body Hamiltonian. In one 

dimension wc have truncated the many-body basis at cut-off energy E„iax = 38 a;. This corresponds to a basis with 
blocks of up to 80524 many-body states. In two dimensions the number of overlap integrals for a given Emax is 
considerably higher than in one dimension. Moreover, the generation of matrix elements turns out to be very time 
consuming. Wc have therefore stopped at Emax = Uuj with a basis with blocks of up to 4532 states in two dimensions. 
Energies are measured in units of the oscillator frequency and lengths in units of the typical oscillator width 1 / ^mio. 
g refers to the corresponding dimensionless interaction strength {giD,92D) as obtained in the appendix. 
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FIG. 1. Convergence of the ground state energy. We plot the ratio of the ground state energy Eq for the given cut-off Emax and the 
energy By""^ found using the highest Emax ■ The upper graph is for one dimension and the lower for two dimensions. 



In Fig. 1 we show the ground state energy as function of the cut-off energy Emax- The plotted quantity is the ratio 
of the ground state energy for the given E^ax and the energy found using the highest E^^x considered. The upper 
graph is for one dimension and the lower for two dimensions. In one dimension, the relative difference in ground state 
energy going from 36 a; to 38 a; is seen to be about 0.05% or less. Going from 12 a; to 14 a; in two dimensions, the shift 
is about 0.2%. The convergence of the excitation energies is found to be considerably faster for the lowest levels, but 
is slower at high energy. This will be seen below in the deviation of the interaction independent energies from the 
exact values mo. 
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FIG. 2. Ground state energy Eo per partiele in one dimension measured in units of w as function of interaetion strength g for particle 
numbers N = 10, 12, 15, 20, 30 and 40 starting from below. The results of the many-body calculation converge towards the mean field 
result as the particle number increases. 



A. Ground state energy 

In this subsection we show the effect of the interaction on the energy of the ground state. In all figures the ideal 
gas ground state energy is set to zero. 

Linear Bose gas. Fig. 2 shows the ground state energy in one dimension, both in the many-body description and 

in the Gross-Pitacvskii approximation. For a given Ng the energy resulting from the many-body calculation is 
seen to converge towards the mean field value as the particle number N is increased. The discrepancy between the 
approximations increases with Ng. 
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FIG. 3. Effective coupling Ng in one (left) and two (right) dimensions. The corresponding values of Ngo are found at 1/N = 0. 
Effective couplings based on the energy and the core density of the ground state are compared. 
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The discrepancies found between the two approximations might to some extent be removed by the introduction of 
an effective couphng constant in the mean field theory. For a set of couphngs go and particle numbers N we find the 
rescaled coupling which in the GP-approximation yields the corresponding energy given by the many-body calculation. 
In Fig. 3 this rescaled coupling multiplied by N is plotted (filled circles) as a function of 1/N for various values of 
the bare coupling. The corresponding value of Ngo is found at 1/iV = 0. For the range of parameters considered, the 
lines look very linear. Below, we compare these results with the effective coupling based on calculations of the core 
density in the ground state. 



Planar Bose gas. The insertion in Fig. 4 shows the Gross-Pitaevskii ground state energy in two dimensions as 
function of Ng. The energy falls rapidly with decreasing, negative coupling. With a basis of up to 28 oscillator 
eigenf unctions, we have not been able to reach convergence near Eq/N = —u. It seems though that dEQ/d{Ng) 
approaches infinity as we get closer to this point, thus indicating that no stable ground state may be formed for 
stronger attractive couplings. This is agreement with Pitaevskii's analysis of the system |24 , where it is shown that 



an attractive interaction leads to collapse of the gas for energies more than uj below the ideal gas ground state energy. 
Similar conclusions are drawn in [E5|. 
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FIG. 4. Ground state energy Eq per particle in two dimensions measured in units of u) as function of interaction strength g for particle 
numbers N = 10, 15 and 20 starting from below. The results of the many-body calculation are compared to the mean field approximation 
(Gross-Pitaevskii), which is approached as the particle number increases. The insertion shows the energy in the mean field approximation 
over a broader range of the interaction strength. 

The Gross-Pitaevskii ground state energy is in Fig. 4 also plotted together with the results of the many-body 
calculation for N — 10, 15 and 20. The mean field approximation is approached as N increases. We have calculated 
the effective coupling Ng which, used in the GP-approximation reproduces the energy obtained in the many-body 
calculation. This effective coupling is plotted as filled circles in the right graph of Fig. 3. The results are all for the 
highest cut-off energy Emax — 14a;. We again find that Ng depends nearly linearly on 1/N . 
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FIG. 5. Excitation energies (E — Eo)/lo in a one-dimensional system of 30 bosons as function of tlie interaction strength g. 

B. Excitations 

In Section ^ we saw that the excitation energies must consist of ladders with energy spacing lo. This is due to the 
existence of the collective raising operator . In two dimensions the breathing mode leads to additional ladders with 
spacing 2 to. The bottom of each of the ladders must however be determined by other means, and here we present the 
results of the numerical diagonalization of the many-body Hamiltonian. 




(E-Eo)/o) 

FIG. 6. Summed density of states in a one-dimensional (right curves, nio) and two-dimensional (left curves, 71213) system of 20 and 
15 bosons respectively. The interaction strength is in the range Ng = 0, 0.2, . . . , 2 in one dimension and Ng = — 1, —0.5, . . . , 2.5 in two 
dimensions. Ng increases from right to left, and negative interaction strengths are indicated by dashed lines. 



Linear Base gas. In Fig. 5 the lowest excitation energies in one dimension are given as function of Ng for the case 
of 30 particles. The degeneracy at zero coupling is split by the interaction. For each level in the ideal case, one level 
remains independent of the coupling. These levels correspond to excitations described by from the ground state. 
The others spread out below, and above each of these we find the mentioned ladder of levels. Similar results are 
obtained for A'' = 10, 12, 15, 20 and 40. The relative shift in the energies going from A'^ = 10 to A'' = 40 with a given 
Ng is less than 1% for the range of parameters presented in Fig. 5. 

The splitting of the degeneracy by the interaction leads to a more even distribution of energy levels than in the free 
case. As an illustration of this. Fig. 6 shows the number of energy levels nmiE) (right curves) below a given energy 
E measured from the ground state. We have taken N = 20, and the coupling ranges from Ng = — 2 in steps of 0.2. 
This summed density of states is indeed smoothed by the presence of the interaction. Note that the curves coincide 
at points {E—Eq)/w = n. This happens because all interaction-dependent excitation energies lie in the gap below the 
corresponding undisturbed level for the considered values of the coupling. 
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FIG. 7. Excitation energies {E — Eq)/uj in a two-dimensional system of 20 bosons as function of the interaction strength. L is the 
angular momentum. 

Planar Bose gas. The lowest energy levels of a planar system with 20 bosons arc given in Fig. 7 as function of the 
dimensionless coupling. For each level with non-zero angular momentum L, there is a level with the same energy and 
angular momentum —L. The degeneracy at <? = is partially lifted for non-zero g, but again there are levels which 
are unaffected by the interaction. These are reached by the two dipole excitations and the breathing mode excitation 
from the ground state. The modes also give rise to a set of ladders above the interaction dependent levels. Since the 
energies of all three interaction independent modes are multiples of o), there is a high degree of degeneracy even for 
5 7^ in this system. 

In Fig. 6 the corresponding summed density of states n2D{E) (left curves) is plotted for 15 bosons with the coupling 
running from Ng = —1 to 2 in steps of 0.5. Levels for both positive and negative angular momentum are included. 
The step-like form of the curve for the free case is again smoothed by the interaction. 

For stronger couplings we do not quite reproduce the independence of the interaction. This is due to the truncation 
of the basis of states. However, since we know that these levels must approach {E—Eq)/u> = integer as the basis set 
is increased, the deviation gives a good indication of the accuracy of the results. In two dimensions we find that the 
convergence is better for those levels involving only dipole excitations than for those which are reached (partly) by 
the breathing mode. As discussed above, this is an expected effect of the existence of the energy cut-off. 
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C. Ground state density 
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FIG. 8. Normalized ground state density p in one dimension. The density is plotted as function of the dimensionless radius x = ry/mU. 
The solid curves give the density in the Gross-Pitaevskii approximation for Ng = —1,0,1,2,5,10 and 50. The dotted curves show the 
result of the Thomas-Fermi approximation for repulsive interaction. The density obtained in the exact many-body calculation with 30 
particles is shown as dashed lines for Ng = 1 and 2. 

The ground state can in the GP-approximation be written as a superposition of oscillator functions with even parity 
or zero angular momentum in one and two dimensions respectively. In the many-body description this is no longer 
true, and single particle states with other quantum numbers will contribute. We here show how this affects the ground 
state density. 

For large gN the interaction energy dominates the kinetic energy, and one may neglect the kinetic term in the 
GP-equation. This is the Thomas- Fermi (TF) approximation 0. The solution is now simply 



1 



2gN 



(Ai-K.(r))0(M-Kx(r)) 



(4.1) 



This describes a condensate density which vanishes outside the radius R — y/2^/niuj^ . We here compare the normal- 
ized density p in the TF-approximation with the exact solution of the Gross-Pitaevskii approximation, and also with 
the full many-body calculation. The TF-approximation in three dimensions has been studied earlier in ||^ . 

In Fig. 8 and Fig. 9 we plot the normalized ground state density p{x) as function of the dimensionless radius 
X = ry^muj in one and two dimensions, respectively. The Thomas-Fermi approximation is seen to reproduce core 
density p{x = 0) for interaction strengths above Ng ~ 5 in one dimension and above Ng ^ 10 in two dimensions. 
The approximation gets worse as the distance from the core is increased, where the density has a sharp cut-off in 
the Thomas-Fermi approximation. For the lowest values of the interaction strength we compare the Gross-Pitaevskii 
approximation to the full many-body calculation with = 30 and = 20 in one and two dimensions, respectively. 
The GP-approximation yields quite good results, especially at large distances, but it tends to overestimate the effect 
of the interaction. Still, the mean field theory seems to work better for the density than for the ground state energy. 
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FIG. 9. Normalized ground state density p in two dimensions. The density is plotted as function of the dimensionless radius x = ry/niuj. 
The solid curves give the density in the Gross-Pitaevskii approximation for Ng = —1,0,1,2,5,10 and 50. The dotted curves show the 
result of the Thomas-Fermi approximation for repulsive interaction Ng = 5, 10 and 50. The density obtained in the exact many-body 
calculation with 20 particles is shown as dashed lines for Ng < 5. 

We may again try to compensate for the deviation by introducing an effective interaction strength, as done for the 
energy calculations. The corresponding values of Ng are shown as triangles in Fig. 3. The one- and two-dimensional 
results are for Emax/i^ — 38 and Emax/^ = 14, respectively. 

The effective couplings are different from those based on the energy calculation, and the relative difference increases 
with NgQ. The agreement is however quite good for Ngo = 0.5 and 1 in two dimensions. But, in general, a simple 
rescaling of the interaction cannot account for all many-body effects in the ground state. 



Having obtained the lowest excitation energies of the system of trapped bosons, it is straightforward to find the 
partition function and the thermodynamical quantities one may derive thereof. The validity of these results is limited 
to low temperatures T Emax- Just how low the temperature must be will be clear in the examples which follow. 

Let us first consider the one-dimensional system. With 20 particles and a cut-off at Emax — 26 a;, wc find the 
specific heat plotted in the left graph of Fig. 10. The coupling ranges from Ng = to 2 in steps of 0.5, starting 
from below at low temperatures. The dotted line is the specific heat for 20 non-interacting bosons calculated in the 
canonical ensemble. We also show the result of a calculation in the grand canonical ensemble. For particle numbers 
as low as this, the two ensembles predict somewhat different specific heats. The exact curve agrees with the curve 
for g = up to a temperature around 1.5 a;. At higher temperatures the many-body calculation breaks down. For 
temperatures where the results are valid we find that the repulsive interaction yields a higher specific heat than in 
the case of non-interacting particles. This happens because the interaction brings most of the energy levels closer to 
the ground state. The effect of the interaction grows with increasing temperature in this regime. 

We have also calculated the specific heat for a two-dimensional system of 15 bosons. The cut-off is here E^ax — 14 a;, 
and the basis thus contains states with angular momentum up to L = ±14. The specific heat is shown in the right 
graph of Fig. 10. The calculation now breaks down around T/lo — 1. For comparison we have included the specific 
heat of 15 non- interacting bosons as given in the grand canonical ensemble. Again a repulsive interaction increases 
the specific heat, while it is decreased in the attractive case. 

There is no Bose-Einstein condensation for an ideal, homogeneo us gas in two dimsneions. In the presence of 
a harmonic trapping potential, condensation may however occure [P7| , p8| . The condensation temperature To = 
hujyj N /C,{2) persists in the thermodynamic limit, which is found by taking A'^ — > oo , a; — > while keeping the average 
particle density p ^ Nuj"^ constant. Recently, MuUin has studied a two-dimensional system with attractive 
interaction in the thermodynamic limit. He finds that the normal state has a transition temperature T < Tq, but not 



D. Specific heat 
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to a Bose condensed phase. This agrees with the earUer findings of Shevchenko ||2^. Keeping, instead N and lu finite, 
one should find a non-zero temperature at which condensation sets in. Although one might study this by the method 
used here, we can not derive any decisive results as far as the condensation temperature is concerned. For this one 
would need excitation energies going considerably higher than what we have considered. 




T/co T/co 

FIG. 10. Specific heat Cv per particle measured in units of kg for a one-dimensional (left) and two-dimensional (right) system of 20 
and 15 bosons respectively. The specific heat is plotted as function of the reduced temperature T/uj. In one dimension the interaction 
strength runs from Ng = to 2 and in two dimensions from -1 to 2.5. Ng increases from below in steps of 0.5. The many-body results 
(mb) are compared to the curve for the ideal gas as given in the grand canonical ensemble (GC). In one dimension we also show the exact 
result for the ideal gas in the canonical ensemble (EC). Notice the breakdown at higher temperatures due to the cutoff in the energy 
spectra. 



V. DISCUSSION AND CONCLUSION 



We have in this paper presented results of calculations in the configuration interaction approximation on systems 
with 10 to 40 interacting bosons confined to harmonic traps in one and two dimensions. By use of a set of many-body 
states and decomposition of the fields in harmonic modes, we have diagonalized the Hamiltonian numerically. The 
resulting ground state properties have been compared to the predictions of the Gross-Pitaevskii equation. The largest 
difference is found for the ground state energy, while for the ground state density, the mean field predictions are 
considerably more accurate, at least for larger Ng. 

We have tried to compensate by introducing an effective interaction strength in the mean-field description. This 
has been done by computing that interaction strength which, when applied in the Gross-Pitaevskii equation yields 
the same ground state energy and core density as found in the full many-body calculation. In two dimensions, and for 
Ngo = 0.5 and 1 the effective couplings obtained from matching energy and from matching core density agree quite 
well. But, in general, rescaling based on ground state energy and core density result in different effective couplings. 
Thus, this simple rescaling cannot account for all many-body effects in the ground state. 

We have also considered the lowest excitation energies. In both one and two dimensions we have reproduced the 
interaction independent levels which must be present in these systems |l^ , [l9|] . The interaction breaks the degeneracy 
in the non-interacting gas completely in one dimension. In two dimensions the two dipole modes and the breathing 
mode assure a high degree of degeneracy even for non-zero interaction. We have calculated the specific heat based 
on these excitation energies. For low temperatures, where the results are valid, the repulsive interaction increases the 
specific heat, while it is decreased in the attractive case. 

It would be interesting to compare the obtained excitation energies with the results of calculations based on the 
Bogoliubov or the Hartree-Fock approximation. One might also use the obtained density of states as a basis for 
calculations in the microcanonical ensemble. Finally the parameter range might be extended by choosing a more 
suitable basis of single particle states, such as the Hartree-Fock states. 
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APPENDIX A: EFFECTIVE COUPLING IN LOWER DIMENSIONS 

Inspired by proposed experiments in lower dimensions |^,^ , we will here consider the effect of having a strongly 
anisotropic trapping potential. Assuming uj^ ^ i^x.y or oj^ 3> ^x,y the system should effectively be one- or two- 
dimensional. When the larger oscillator energy dominates all other energies in the system, we may assume that there 
will be no excitations in this direction. This means that the field operator may be written as 

^{^iV.z) = ^(/)o(a;)0o(y)0fc(2;)aoofc (Al) 
fc 

for the one dimensional case, and 

^(a;,y, z) = ^(/)o(z)(/)fc(a;,y)aofe (A2) 
fc 

for the two-dimensional case. We may then integrate over the spatial varia bles corresponding to the stronger oscillation 



frequencies. This results in a rescaling of the coupling constant g in Eq. (2.2) 



92D=93Dsj^^ and .gu = .gsD ^ ■ (A3) 

Measuring energy in units of the remaining oscillator frequency lo and lengths in units of the oscillator length 1/ ^rtiLj 
we obtain the dimensionless coupling constants 



ff2D = V27r— and 5in = tt,/ — 2- , (A4) 

Qz V '^z ax 



ay 

where = ^/h/muTi and a = mg^D /'^ir'h? ■ We have here reintroduced fi. The effective coupling is seen to increase 
with the strength of the stronger oscillator frequency. This is an effect of reducing the length in which the particle 
density may extend in the dimensions integrated out. 

As an example of a possible physical realization, we mention that Ketterle and van Druten have suggested an 
experiment with Lj^l^xy ~ 10~^ and uj^y ^ 2{){)nK . With a scattering length of typically 100 Bohr radii, this yields 
an effective one-dimensional coupling ^ 1. This corresponds to a strong coupling system for any finite number of 
particles. Other configurations may yield a weaker effective coupling. 
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